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We use some fractal analysis methods to study river flow fluctuations. The result of the Multi- 
fractal Detrended Fluctuation Analysis (MF-DFA) shows that there are two crossover timescales 
at six ~ 12 and S2 X ~ 130 months in the fluctuation function. We discuss how the existence of 
the crossover timescales are related to a sinusoidal trend. The first crossover is due to the seasonal 
trend and the value of second ones is approximately equal to the well known cycle of sun activity. 
Using Fourier detrended fluctuation analysis, the sinusoidal trend is eliminated. The value of Hurst 
exponent of the runoff water of rivers without the sinusoidal trend shows a long range correlation 
behavior. For the Daugava river the value of Hurst exponent is 0.52 ± 0.01 and also we find that 
these fluctuations have multifractal nature. Comparing the MF-DFA results for the remaining data 
set of Daugava river to those for shuffled and surrogate series, we conclude that its multifractal 
nature is almost entirely due to the broadness of probability density function. 
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I. INTRODUCTION 

Interpretation and estimation of climate change has 
been one of the main research areas in science [1-3] . The 
climate system often exhibits irregular and complex be- 
havior. Although the climate system is driven by the 
well-defined seasonal periodicity, it is also a subject to 
unpredictable perturbations which can lead to extreme 
climate events. Indeed the climate is a dynamical system 
influenced by immense factors, such as solar radiation or 
the topography of the surface of the solid earth, etc. All 
factors that control the trajectory of climate have enor- 
mously large phase space, thus we have to analysis it 
with stochastic tools. Several recent statistical studies 
have shown that a remarkably wide variety of natural 
systems display fluctuations that may be characterized 
by long-range power-law correlations. Such correlations 
hint toward fractal geometry of the underlying dynamical 
system. Existence and determination of power-law cor- 
relations would help to quantify the underlying process 
dynamics [4,5]. 

The analysis of river flows has a long history, neverthe- 
less some important issues have been lost. Here, we study 
one component of the climate system, the river flux, by 
using the novel approach in the fractal analysis like De- 
trended Fluctuation Analysis, Fourier-Detrended Fluctu- 
ation Analysis and Scaled Windowed Variance Analysis 
Methods. The statistical and fractal analysis of river 
flows should be an important issue in the geophysics and 
hydrological systems to recognize the influence of envi- 
ronmental conditions and to detect relative effects. A 
set of most important results which can be given by us- 
ing statistical tools is as follows: a concept of scale self- 
similarity for the topography of Earth's surface [6], the 



hydraulic-geometric similarity of river system and floods 
forced by the heavy rain [7,8], etc. Already more than 
half a century ago the engineer Hurst found that runoff 
records from various rivers exhibit 'long-range statisti- 
cal dependencies' [9]. Later, such long-term correlated 
fluctuation behavior has also been reported for many 
other geophysical records including precipitation data 
[6,10,11]. These original approaches exclusively focused 
on the absolute values or the variances of the full dis- 
tribution of the fluctuations, which can be regarded as 
the first and second moments of detrended fluctuation 
analysis [6,9,10,12]. In the last decade it has been real- 
ized that a multifractal description is required for a full 
characterization of the runoff records [13,14]. This mul- 
tifractal description of the records can be regarded as a 
'fingerprint' for each station or river, which, among other 
things, can serve as an efficient non-trivial test bed for 
the state-of-the-art precipitation-runoff models. 

River flow can be characterized by several general fea- 
tures. As a result of the periodicity in precipitation, river 
flow has also strong seasonal periodicity. The seasonal 
cycle of river flow is asymmetric; i.e., river flow increases 
rapidly (usually during late winter and spring) and de- 
creases gradually (toward the end of the autumn). The 
fluctuations in river flow arc large for large river flow and 
small for small river flow [5] . It is important to note that 
unlike other climate components, river flow may has a 
direct impact of human activity, like damming, use of 
river water for agriculture, etc., a fact which makes the 
river flow data more difficult to study. The fluctuations 
in river flow are of special interest since they are directly 
linked to floods and droughts. There are several inter- 
esting characteristics of river flow fluctuations: (i) the 
river flow fluctuations have power law tails in the proba- 



bility distribution [15,16], (ii) the river flow fluctuations 
are long-range correlated [9,17,18], and (iii) river flow 
fluctuations are multifractal [13,14,19]. The scaling laws 
may improve the statistical prediction of extreme changes 
in river flow [20]. Recently connection between volatil- 
ity and nonlincarity has originally been established by 
[21,22], the degree of non-linearity has been checked us- 
ing the volatility series, also a simple model of river fluc- 
tuations has been determined [4,5]. More recently the 
annual runoff for the Ukrainian and Moldavia's rivers 
and reveal scale invariance for distribution of this vari- 
able have been investigated by using statistical parame- 
ters such as arithmetic average, coefficients of variation, 
skewness, and auto-correlation [23]. In all of the pre- 
vious researches, the contribution sinusoidal trends on 
the creation of crossover in the results of fractal analysis 
and the multifractal nature have been lost. The effect of 
nonstationarity on the detrended fluctuation analysis has 
been investigated [24] . In addition, the effects of periodic 
trends on the fractal scaling properties of a time series 
have been investigated moreover in some paper a relation 
between amplitude and the period of the periodic trends 
and the existence of crossover in the Detrended Fluctu- 
ations scaling function have been demonstrated [24]. So 
the main purpose of this paper are the investigation of 
the effects of seasonal trend on the multifractal analysis 
of flow fluctuations and determination of the source of 
multifractality in data . For completeness of this inves- 
tigation and to get the deep insight of the contribution 
of sun activity in the statistical properties of river flow, 
we compare the recent fractal analysis results of sunspots 
[25] with current analysis. The Sunspot number was col- 
lected by the Sunspot Index Data Center (SIDC) [26]. 
Due to the stochastic nature of river flow, it is probable 
that the sun activity may affects on the duct of river, 
so to demonstrate the presence of any correlations we 
compare results of river flow and sun activity extracted 
by various fractal analysis. It was well-known that the 
statistical properties of every rivers depend on very im- 
portant reasons which affect on flow fluctuations, so one 
can not expect that the sun activity has a same reason- 
able effect on different rivers. 

In addition we would like to characterize the complex 
behavior of the monthly runoff for the Daugava river fluc- 
tuations through the computation of the signal parame- 
ters - scaling exponents - which quantifies the correlation 
exponents and multifractality of the signal. We inves- 
tigate the correlation behavior of duct river time series 
which is governed by power-law. 

The original Daugava river data source is Latvian 
Environmental Geological and Meteorological Agency 
database. They describe water flow through hydroelec- 
tric power station near Kegums, Latvia. Dimension of 
the data is m 3 /s. Other data which are used here are 
from National Water Information System: Web Inter- 
face [27]. As shown in the upper panel of Figure 1, the 



duct water of Daugava river scries has a sinusoidal trend, 
with a dominant frequency. These trends should involve 
the seasonal and other physical reasons in natural phe- 
nomenon. Because of the complexity nature of river flow 
series, and due to the finitcness of the available data sam- 
ple, we should apply some methods which are insensitive 
to non-stationarities like trends. To eliminate the effect 
of sinusoidal trend, we apply the Fourier Detrended Fluc- 
tuation Analysis (F-DFA) [28,29]. After elimination of 
the trend we use the Multifractal Detrended Fluctua- 
tion Analysis (MF-DFA) to analysis the data set. The 
MF-DFA methods are the modified version of detrended 
fluctuation analysis (DFA) to detect multifractal prop- 
erties of time series. The detrended fluctuation analysis 
(DFA) method introduced by Peng et al. [30] has be- 
came a widely-used technique for the determination of 
(multi-) fractal scaling properties and the detection of 
long-range correlations in noisy, nonstationary time se- 
ries [24,30-33]. It has successfully been applied to di- 
verse fields such as DNA sequences [30,34], heart rate 
dynamics [35-37], neuron spiking [38], human gait [39], 
long-time weather records [40], cloud structure [41], ge- 
ology [42], ethnology [43], economical time series [44,45], 
solid state physics [46], sunspot time series [25] and cos- 
mic microwave background radiation [47] . 

This paper is organized as follows: In Section II we 
describe the MF-DFA, F-DFA and Scale Windowed Vari- 
ance (SWV) methods in detail and show that the scal- 
ing exponents determined via the MF-DFA method are 
identical to those obtained by the standard multifractal 
formalism based on partition functions. We eliminate the 
sinusoidal trend via the F-DFA technique in Section III 
and investigate the multifractal nature of the remaining 
fluctuation, we use certain fractal analysis approaches 
such as, the Multifractal Detrended Fluctuation Analy- 
sis (MF-DFA), and Scaled Windowed Variance (SWV) 
to analysis the data set, The DFA result of sun activity 
and river flows arc compared together. In Section IV, 
we examine the source of multifractality in duct water of 
Daugava river data by comparison the MF-DFA results 
for remaining data set to those obtained via the MF-DFA 
for shuffled and surrogate series. Section V closes with a 
discussion of the present results. 



II. FRACTAL ANALYSIS METHODS 

In this section we introduce three methods to investi- 
gate the fractal properties of stochastic processes. 

A. Multifractal Detrended Fluctuation Analysis 

The simplest type of the multifractal analysis is based 
upon the standard partition function multifractal for- 
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FIG. 1. Observed flux series of Daugava river (upper panel) 
and Sunspot number (lower panel) as a function of time. 



malism, which has been developed for the multifractal 
characterization of normalized, stationary measurements 
[48-51]. Unfortunately, this standard formalism does not 
give correct results for nonstationary time series that are 
affected by trends or that cannot be normalized. Thus, in 
the early 1990s an improved multifractal formalism has 
been developed, the wavelet transform modulus maxima 
(WTMM) method [52], which is based on the wavelet 
analysis and involves tracing the maxima lines in the 
continuous wavelet transform over all scales. The other 
method, the multifractal detrcndcd fluctuation analysis 
(MF-DFA) , is based on the identification of scaling of the 
qth-order moments depending on the signal length and is 
generalization of the standard DFA using only the second 
moment q = 2. 

The MF-DFA does not require the modulus maxima 
procedure in contrast WTMM method, and hence does 
not require more effort in programming and computing 
than the conventional DFA. On the other hand, often 
experimental data are affected by non-stationaritics like 



trends, which have to be well distinguished from the in- 
trinsic fluctuations of the system in order to find the cor- 
rect scaling behavior of the fluctuations. In addition very 
often we do not know the reasons for underlying trends 
in collected data and even worse we do not know the 
scales of the underlying trends, also, usually the avail- 
able record data is small. For the reliable detection of 
correlations, it is essential to distinguish trends from the 
fluctuations intrinsic in the data. Hurst rescaled-rangc 
analysis [10] and other non-detrending methods work well 
if the records are long and do not involve trends. But if 
trends are present in the data, they might give wrong 
results. Detrcndcd fluctuation analysis (DFA) is a well- 
established method for determining the scaling behavior 
of noisy data in the presence of trends without knowing 
their origin and shape [30,36,53-55]. Also DFA scaling 
results are not immune to different trends and to differ- 
ent artifacts such as spikes, missing segments of data etc. 
[33,56,57]. 

In spite of many abilities of this method, in some 
cases is encountered with problem and gives wrong re- 
sults. DFA method can only determine positive Hurst 
exponent, H, and gives an inaccurate results for strongly 
anti-correlated record data when H is close to zero. To 
avoid this situation, in such case one should use the inte- 
grated data. This signal is so-called double profiled data 
set. The corresponding Hurst exponent using this way is 
H = H — 1, here H is derived from DFA method for the 
double profiled signals [25,58]. According to the recent 
exploration in ref. [32], a deviation in the DFA results 
occurs in very short records and in the small regime of 
data in each window mentioned in the forthcoming sub- 
section. The modified version of MF-DFA will be used 
for these cases [32]. In this paper also we will introduce 
this method and apply to infer correct exponent for river 
flow fluctuations. 



B. Description of the MF-DFA 

The modified multifractal DFA (MF-DFA) procedure 
consists of five steps. The first three steps are essentially 
identical to the conventional DFA procedure (see e. g. 
[24,30-33]). Suppose that Xk is a series of length A, and 
that this series is of compact support, i.e. Xk = for an 
insignificant fraction of the values only. 
• Step 1: Determine the "profile" 



Y(i) EE J2 [Xk - (X)} , 



1,...,N. 



(1) 



fc=l 



Subtraction of the mean (x) is not compulsory, since it 
would be eliminated by the later detrending in the third 
step. 

• Step 2: Divide the profile Y(i) into N s = int (iV/s) non- 
overlapping segments of equal lengths s. Since the length 



TV of the series is often not a multiple of the considered 
time scale s, a short part at the end of the profile may 
remain. In order not to disregard this part of the series, 
the same procedure is repeated starting from the opposite 
end. Thereby, 2N S segments are obtained altogether. 
• Step 3: Calculate the local trend for each of the 2N S 
segments by a least-square fit of the series. Then deter- 
mine the variance 



F 2 ( S ,v)^-J2{Y[(u-l)s + i}-y v {i)} 2 



(2) 



for each segment u, v = 1, . . . , N a and 

F 2 (s, v) = - V {Y[N -{v- N s )s + i] - y^)} 2 , (3) 
s 

i=l 

for v = N s + 1, . . . , 2N S . Here, y v (i) is the fitting poly- 
nomial in segment v. Linear, quadratic, cubic, or higher 
order polynomials can be used in the fitting procedure 
(conventionally called DFAl, DFA2, DFA3, . . .) [30,37]. 
Since the detrending of the time series is done by the 
subtraction of the polynomial fits from the profile, dif- 
ferent order DFA differ in their capability of eliminating 
trends in the series. In (MF-)DFAm [mth order (MF- 
)DFA] trends of order m in the profile (or, equivalcntly, 
of order m— 1 in the original series) are eliminated. Thus 
a comparison of the results for different orders of DFA al- 
lows one to estimate the type of the polynomial trend in 
the time series [24,32]. 

• Step 4'- Average over all segments to obtain the g-th 
order fluctuation function, defined as: 



F q {s) 



1 



2N a 
u—1 



1/9 



q/2 



(4) 



where, in general, the index variable q can take any real 
value except zero. For q = 2, the standard DFA proce- 
dure is retrieved. Generally we are interested in how the 
generalized q dependent fluctuation functions F q (s) de- 
pend on the time scale s for different values of q. Hence, 
we must repeat steps 2, 3 and 4 for several time scales 
s. It is apparent that F q (s) will increase with increasing 
s. Of course, F q (s) depends on the DFA order m. By 
construction, F q {s) is only defined for ,s > m + 2. 
• Step 5: Determine the scaling behavior of the fluctua- 
tion functions by analyzing log- log plots of F q (s) versus s 
for each value of q. If the series Xi are long-range power- 
law correlated, F q {s) increases, for large values of s, as a 
power-law, 



F q (s) 



Mi) 



(5) 



In general, the exponent h(q) may depend on q. For 
stationary time series such as fGn (fractional Gaussian 
noise), Y(i) in Eq. 1, will be a fBm (fractional Brownian 



motion) signal, so, < h(q = 2) < 1.0 (see the appendix 
for more details). The exponent h(2) is identical to the 
well-known Hurst exponent H [30,31,48]. Also for a non- 
stationary signal, such as fBm noise, Y(i) in Eq. 1, will 
be a sum of fBm signal, so the corresponding scaling ex- 
ponent of F q (s) is identified by h(q = 2) > 1.0 [30,59]. 
In this case the relation between the exponents h(2) and 
H will be H = h(q = 2) — 1 (see appendix of [25]). 
The exponent h(q) is known as generalized Hurst expo- 
nent. The auto-correlation function can be characterized 
by a power law C(s) = (rifc^fc+s) ~ s -7 with exponent 
7 = 2 — 2H. Its power spectra can be characterized by 
S(ui) ~ lq~P with frequency ui and [3 = 2H — 1. In the 
nonstationary case, correlation function is (see appendix 
for more details): 

\2H 



C(iJ) 



iUj) 



J\ 



(0) 



and corresponding power spectrum scaling is (3 = 2H + 1 
[25,30,59]. 

For monofractal time series, h(q) is independent of q, 
since the scaling behavior of the variances F 2 (s,v) is 
identical for all segments v, and the averaging procedure 
in Eq. (4) will just give this identical scaling behavior 
for all values of q. If we consider positive values of q, 
the segments v with large variance F 2 (s,v) (i. e. large 
deviations from the corresponding fit) will dominate the 
average F g (s). Thus, for positive values of q, h{q) de- 
scribes the scaling behavior of the segments with large 
fluctuations. For negative values of q, the segments v 
with small variance F 2 (s,v) will dominate the average 
F q (s). Hence, for negative values of q, h(q) describes the 
scaling behavior of the segments with small fluctuations 
[60]. 



1. Relation to standard multifractal analysis 

For a stationary, normalized scries the multifractal 
scaling exponents h(q) defined in Eq. (5) are directly re- 
lated to the scaling exponents r(q) defined by the stan- 
dard partition function-based multifractal formalism as 
shown below. Suppose that the series Xk of length N is 
a stationary, normalized sequence. Then the detrending 
procedure in step 3 of the MF-DFA method is not re- 
quired, since no trend has to be eliminated. Thus, the 
DFA can be replaced by the standard Fluctuation Anal- 
ysis (FA), which is identical to the DFA except for a 
simplified definition of the variance for each segment u, 
v = 1, . . . , N s . Step 3 now becomes [see Eq. (2)]: 



F 2 A (s,v) = [Y(vs)-Y((v-l)s)Y 



(7) 

Inserting this simplified definition into Eq. (4) and using 
Eq. (5), we obtain 

1/9 



2JV S 

— y 

2N S ^ 



Y(ys) -Y((u- l)s)\ q 



Mi) 



(8) 



For simplicity we can assume that the length TV of the 
series is an integer multiple of the scale s, obtaining N s = 
N/s and therefore 



N/s 



Y,\Yiy*)-Y{iy-i)*)\ q 



v qKq)-i 



(9) 



This corresponds to the multifractal formalism used e. g. 
in [49,51]. In fact, a hierarchy of exponents H q simi- 
lar to our h(q) has been introduced based on Eq. (9) in 
[49]. In order to relate also to the standard textbook 
box counting formalism [48,50], we employ the defini- 
tion of the profile in Eq. (1). It is evident that the term 
Y{vs) — Y((v — l)s) in Eq. (9) is identical to the sum 
of the numbers Xk within each segment v of size s. This 
sum is known as the box probability p s (y) in the standard 
multifractal formalism for normalized series Xk, 

us 

Ps {v)= Yl x k =Y(vs)-Y((v-l)s). (10) 

fe=(i/-l)s+l 

The scaling exponent r(q) is usually defined via the par- 
tition function Z q (s), 



Z q {s) 



N/s 

£ 

i/=i 



Pa{v)\ 



e r(9) 



(11) 



where q is a real parameter as in the MF-DFA method, 
discussed above. Using Eq. (10) we see that Eq. (11) is 
identical to Eq. (9), and obtain analytically the relation 
between the two sets of multifractal scaling exponents, 

T(q) = qh(q) - 1. (12) 

Thus, we observe that h(q) defined in Eq. (5) for the 
MF-DFA is directly related to the classical multifractal 
scaling exponents r(q). Note that h{q) is different from 
the generalized multifractal dimensions 

r{q) qh(q) 



D{q) 



1 



(13) 
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that are used instead of r(q) in some papers. While h(q) 
is independent of q for a monofractal time series, D(q) 
depends on q in this case. Another way to characterize a 
multifractal series is the singularity spectrum f(a), that 
is related to r(q) via a Legendre transform [48,50], 

a = r'(q) and f(a) = qa — T(q). (14) 

Here, a is the singularity strength or Holder exponent, 
while f(a) denotes the dimension of the subset of the 
series that is characterized by a. Using Eq. (12), we can 
directly relate a and /(a) to h(q), 

a = h(q) + qh'(q) and /(a) = q[a - h(q)] + 1. (15) 

A Holder exponent denotes monofractality, while in 
the multifractal case, the different parts of the structure 
are characterized by different values of a, leading to the 
existence of the spectrum f(a). 



C. Fourier-Detrended Fluctuation Analysis 

Many signals in the nature do not distinguished as 
monofractal scaling behavior. In some cases, there ex- 
ist one or more crossover (time) scales s x segregating 
regimes with different scaling exponents e.g long range 
correlation for s<s x and an other type of correlation or 
uncorrelated behavior for s ^> s x [24,32]. In other cases 
investigation of the scaling behavior is more complicated. 
In the presence of different behavior of various moments 
in the MF-DFA method, different scaling exponents are 
required for different parts of the series [33] . Therefore 
one needs a multitude of scaling exponents (multifrac- 
tality) for a full description of the scaling behavior. A 
crossover usually can arise from a change in the corre- 
lation properties of the signal at different time or space 
scales, or can often arise from trends in the data. To 
remove the crossover due to a trend such as sinusoidal 
trends, Fourier-Detrended Fluctuation Analysis (F-DFA) 
is applied. The F-DFA is a modified approach for the 
analysis of low frequency trends added to a noise in time 
series [28,29,61,62]. 

In order to investigate how we can remove trends hav- 
ing a low frequency periodic behavior, we transform data 
record to Fourier space, then we truncate the first few 
coefficient of the Fourier expansion and inverse Fourier 
transform the series. After removing the sinusoidal 
trends we can obtain the fluctuation exponent by using 
the direct calculation of the MF-DFA. If truncation num- 
bers are sufficient, The crossover due to a sinusoidal trend 
in the log- log plot of F q (s) versus s disappears. 



D. Scaled Windowed Variance Analysis 

The Scaled Windowed Variance analysis was developed 
by Cannon et al. (1997) [59]. The profile of data, Y(i), 
is divided into N s = int(iV/s) non-overlapping segments 
of equal lengths s. Then the standard deviation is calcu- 
lated within each interval using the following relation 



sw( s ) = (i£[Y(i)-(Y( s ))] 5 



1/2 



(16) 



The average standard deviation of all windows of length 
s is computed. This computation is repeated over all 
possible interval lengths. The scaled windowed variance 
is related to s by a power law 



swv 



(17) 



III. ANALYSIS OF DATA 

As mentioned in section II, a spurious of correlations 
may be detected if time series is nonstationarity, so di- 



rect calculation of correlation behavior, spectral density 
exponent, fractal dimensions etc., don't give the reliable 
results. It can be checked that the runoff for Daugava 
river is or not nonstationary. One can verified the non- 
stationarity property experimentally by measuring the 
stability of the average and variance in a moving window 
for example with scale s. Figure 2 shows the standard 
deviation of Daugava flow signal versus scale s, is satu- 
rated. Let us determine that whether the data set has 
a sinusoidal trend or not. According to the MF-DFA1 
method. Generalized Hurst exponents h(q) in Eq. (5) can 
be found by analyzing log-log plots of F q (s) versus s for 
each q. Our investigation shows that there are at least 
two crossover time scales s x in the log-log plots of F q (s) 
versus s for every g's. These two crossovers divide F q (s) 
into three regions, as shown in Figure 3 ( for instance 
we took q = 2). Figure 4 shows these crossover also ex- 
ist in the fluctuation function for various moments. The 
existence of these regions is due to the competition be- 
tween noise and sinusoidal trend [24]. For s < s lx , the 
noise has the dominating effect . For s lx < s < s 2x the 
sinusoidal (such as seasonal) trend dominates [24]. The 
values of s± x and S2 X are approximately equal to 12 and 
130 months, respectively The first crossover is clearly 
related to seasonal trend and the second ones is approxi- 
mately equal to the well known cycle of sun activity. This 
shows that in addition of seasonal effect on the river flow, 
the sun activity strictly affects on the feature and fractal 
properties of river flow fluctuations. As shown in Figure 
3, by comparing curves from some rivers and Sunspot we 
can see a certain symmetry in the form of curves. Points 
of inflection are placed closely, but angle of bendings are 
opposite (e.g. for Daugava river). This symmetry indi- 
cates a connection between the activity of sun and flow 
of water in the rivers. In both cases sinusoidal tendency 
have been found in the average part of the curves. Ap- 
parently the sun activity governs the rivers. 

As mentioned before, for very small scales s < si x the 
effect of the sinusoidal trend is not pronounced, indicat- 
ing that in this scale region the signal can be considered 
as noise fluctuating around a constant which is filtered 
out by the MF-DFA1 procedure. In this region the gen- 
eralized DFA1 exponent for used rivers in this paper are 
listed in Figure 5 , where confirms that the process is a 
stationary process with long-range correlation behavior. 

To cancel the sinusoidal trend in MF-DFAI, we apply 
F-DFA method on the present data. We truncate some of 
the first coefficients of the Fourier expansion of the river 
flow fluctuations. According to Figure 6, for eliminating 
the crossover scales, we need to remove approximately 
the first 200 terms of the Fourier expansion. Then, by 
inverse Fourier Transformation, the noise without sinu- 
soidal trend is extracted (see Figure 6) [28,29,61,62]. 

The MF-DFAI results of the remanning new signal just 
for Daugava river are shown in Figure 7. All exponent in 
this figure are driven at as maximum scale as possible in 
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FIG. 2. Behavior of standard deviation of duct of water 
as a function of time scale. It shows that this time series is 
stationary. 



the log-log plots of F q (s) as function of s (sec Appendix 
for more details). The duct water of Daugava river se- 
ries is a multifractal process as indicated by the strong 
q dependence of generalized Hurst exponents and r(q) 
[58] . The q- dependence of the classical multifractal scal- 
ing exponent r(g) has different behaviors for q < and 
q > 0. For positive and negative values of q, the slopes 
of r{q) are 0.45 ± 0.02 and 1.17 ± 0.02, respectively. Ac- 
cording to the value of Hurst exponent, H = 0.52 ± 0.01, 
we find that this series has approximately random be- 
havior. This result is equal to value of Hurst exponent in 
small scale of MF-DFAl of noise with sinusoidal trend. 
The values of derived quantities from MF-DFAI method, 
are given in Table I. Using the SWV we also analysis the 
truncated series, Our result shows that the value of Hurst 
exponent is H = 0.50 ± 0.01, which is in agreement with 
the previous result. 

Usually, in the MF-DFA method, deviation from a 
straight line in the log-log plot of Eq. (5) occurs for small 
scales s. This deviation limits the capability of DFA to 
determine the correct correlation behavior for very short 
scales and in the regime of small s. The modified MF- 
DFA is defined as follows [32]: 



nmoc 



F q (s) 
K q (sY 

F q (s) 



[j^huf ( s ')]2)l/2 s l/2 



([F shuf ( S )] 2 ) 1 /2 s 'l/2 



(for** > 1), 



(18) 



where ([F g shuf (s)] 2 ) 1 / 2 denotes the usual MF-DFA fluctu- 
ation function [defined in Eq. (4)] averaged over several 
configurations of shuffled data taken from the original 
time series, and s' k iV/40. The values of the Hurst 
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FIG. 3. Crossover behavior of log-log plot F (s) versus s 
for some river flows and Sunspot series for q = 2.0. For rivers 
there are at least two crossover time scales in plot of F(s), at 
scales six and S2 X • 
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FIG. 4. Crossover behavior of log-log plot F(s) versus s 
for Daugava river flow for various moments indicated on the 
figure. For river there are at least two crossover time scales in 
plot of F(s), at scales six an d S2x ■ Vertical axis has arbitrary 
unite. 



exponent obtained by modified MF-DFAl methods for 
duct water of Daugava river series is 0.54 ±0.02. The rel- 
ative deviation of the Hurst exponent which is obtained 
by modified MF-DFAl in comparison to MF-DFAl for 
original data is approximately 3.84%. Now the value of 
Hurst exponent ensure us that the runoff fluctuations are 
long-range correlation, so by ignoring the seasonal trend, 
these processes have almost memory. This means that 
the status of runoff water statistically has long memory. 



IV. COMPARISON OF THE MULTIFRACTALITY 
FOR ORIGINAL, SHUFFLED AND SURROGATE 
SERIES 




FIG. 5. The values of Hurst exponent for some famous 



As discussed in the section III the remanning data set 
after the elimination of sinusoidal trend has the multi- 
fractal nature. In this section we are interested in to 
determine the source of multifractality. In general, two 
different types of multifractality in time series can be dis- 
tinguished: (i) Multifractality due to a fatness of proba- 
bility density function (PDF) of the time series. In this 
case the multifractality cannot be removed by shuffling 
the series, (ii) Multifractality due to different correla- 
tions in small and large scale fluctuations. In this case 
the data may have a PDF with finite moments, e. g. 
a Gaussian distribution. Thus the corresponding shuf- 
fled time series will exhibit mono-fractal scaling, since 
all long-range correlations are destroyed by the shuffling 
procedure. If both kinds of multifractality are present, 
the shuffled series will show weaker multifractality than 
the original series. The easiest way to clarify the type of 
multifractality, is by analyzing the corresponding shuf- 
fled and surrogate time series. The shuffling of time se- 



ries destroys the long range correlation, Therefore if the 
multifractality only belongs to the long range correlation, 
we should find h s hui {o) = 0-5- The multifractality nature 
due to the fatness of the PDF signals is not affected by the 
shuffling procedure. On the other hand, to determine the 
multifractality due to the broadness of PDF, the phase of 
discrete fourier transform (DFT) coefficients of the duct 
water of Daugava river time series are replaced with a 
set of pseudo independent distributed uniform (— ir, n) 
quantities in the surrogate method. The correlations in 
the surrogate series do not change, but the probability 
function changes to the Gaussian distribution. If multi- 
fractality in the time series is due to a broad PDF, h(q) 
obtained by the surrogate method will be independent 
of q. If both kinds of multifractality are present in the 
duct water of Daugava river time series, the shuffled and 
surrogate series will show weaker multifractality than the 
original one. 
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FIG. 6. The MF-DFA1 functions F 2 (s) versus the time 
scale s in log-log plot. Original time series m = 0, trunca- 
tion of m = 50, m = 100 and m = 200 first terms. 




FIG. 7. The q dependence of the generalized Hurst expo- 
nent h(q), the corresponding r(q), generalized multifractal di- 
mension D(q) and singularity spectrum f(a) are shown in the 
upper to lower panel respectively for duct water of Daugava 
river series without sinusoidal trend. 



To check the nature of multifractality, we compare the 
fluctuation function F q (s), for the original series ( af- 
ter cancelation of sinusoidal trend) with the result of 
the corresponding shuffled, F| huf (s) and surrogate series 
F| ur (s). Differences between these two fluctuation func- 
tions with the original one, directly indicate the presence 
of long range correlations or broadness of probability den- 
sity function in the original series. These differences can 
be observed in a plot of the ratio F q (s)/F° hut (s) and 
F q (s) I F| ur (s) versus s [58]. Since the anomalous scal- 
ing due to a broad probability density affects F q (s) and 
F| huf (s) in the same way, only multifractality due to cor- 
relations will be observed in F g (s)/F| huf (s). The scaling 
behavior of these ratios are 



F^/F^is) 



s h(q)-h shat (q) _ ghcoAl) 



(19) 



F q (s)/F^(s) ~ s^-W?) = s ''pdf(9)_ (20) 

If only fatness of the PDF is responsible for the multifrac- 
tality, one should find h(q) = /i s huf(?) and h cor (q) = 0. 
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FIG. 8. Generalized Hurst exponent, h(q) as a function of 
q for original, surrogate and shuffled data. 



On the other hand, deviations from h COI (q) = indicates 
the presence of correlations, and q dependence of h COT (q) 
indicates that multifractality is due to the long rage cor- 
relation. If only correlation multifractality is present, one 
finds /i s huf {q) = 0.5. If both distribution and correlation 
multifractality are present, both, h s huf (<?) and h SUI (q) will 
depend on q. The q dependence of the exponent h(q) for 
original, surrogate and shuffled time series are shown in 
Figure 8. The MF-DFAl results of the surrogate and 
shuffled signal are shown in Figures 9 and 10, respec- 
tively. The q dependence of h cm - and fopDF shows that 
the multifractality nature of the duct water of Daugava 
river time series is due to both broadness of the PDF and 
long range correlation. The absolute value of h COI (q) is 
greater than /jpdf(<z); so the multifractality due to the 
fatness is weaker than the multifractality due to the cor- 
relation. The deviation of h sur (q) and /i s huf (?) from h(q) 
can be determined by using % 2 test as follows: 



JV 

E 



[MgO - Mg.)]' 

^ cr((? ? :) 2 + a (q t ) 2 



(21) 



the symbol " o " can be replaced by " sur" and " shuf ' , 
to determine the confidence level of h sur and /i s huf to 
generalized Hurst exponents of original series, respec- 

2 

tively. The value of reduced chi-square xto = (-M" 
is the number of degree of freedom) for shuffled and sur- 
rogate time series are 6.98, 29.97, respectively. On the 
other hand the width of singularity spectrum, f(ot), i.e. 
Aa = a(<7 m i n ) — a((?max) for original, surrogate and shuf- 
fled time series are approximately, 1.086, 0.390 and 1.210 
respectively. These values also show that the multifrac- 
tality due to the broadness of the probability density 
function is dominant [63]. 
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FIG. 9. The q dependence of the generalized Hurst expo- 
nent h(q), the corresponding r(q), generalized multifractal 
dimension D(q) and singularity spectrum f(a) are shown in 
the upper to lower panel respectively for surrogate duct water 
of Daugava river series without sinusoidal trend. 



The values of the generalized Hurst exponent h(q = 
2.0), multifractal scaling r(q = 2) and other related scal- 
ing exponents are indicated in Table I for the original, 
shuffled of duct water of Daugava river series obtained 
with MF-DFAl method are reported in Table I. The val- 
ues of the Hurst exponent obtained by MF-DFAl and 
modified MF-DFAl methods for original, surrogate and 
shuffled duct water of Daugava river series are given in 
Table II. 



V. CONCLUSION 

The MF-DFA method allows us to determine the mul- 
tifractal characterization of the nonstationary and sta- 
tionary time series. We apply the recent method to in- 
vestigate the existence of crossover on the result of MF- 
DFA of river flow fluctuations. The concept of MF-DFA 
of runoff water of rivers can be used to gain deeper insight 
in to the processes occurring in climate and hydrological 




-1 - • 



-2-10 1 2 

q 




0.7 - 

• ; 

0.4 0.5 0.6 0.7 0.8 0.9 1 U 1.2 
a 

FIG. 10. The q dependence of the generalized Hurst ex- 
ponent h(q), the corresponding r(q), generalized multifractal 
dimension D(q) and singularity spectrum /(a) are shown in 
the upper to lower panel respectively for shuffled duct water 
of Daugava river series without sinusoidal trend. 



TABLE I. The values of the Hurst {H), power spectrum 
scaling (/3), auto-correlation scaling (7) exponents and multi- 
fractal scaling for original, surrogate and shuffled of monthly 
duct water of Daugava river series obtained by MF-DFA1. 



Data 


H 


P 


7 


r 


Original 


0.52 ±0.02 


0.04 ± 0.02 


0.96 ±0.02 


0.04 ± 0.02 


Surrogate 


0.51 ±0.02 


0.02 ±0.02 


0.98 ± 0.02 


0.02 ±0.02 


Shuffled 


0.50 ±0.01 


0.00 


1.00 ±0.02 


0.00 



TABLE II. The value of the Hurst exponent using 
MF-DFA1 and modified MF-DFA1 for the original, shuffled 
and surrogate of monthly duct water of Daugava river series. 



Method 


Original data 


Surrogate 


Shuffled 


MF-DFA1 


0.52 ±0.01 


0.51 ±0.01 


0.50 ±0.01 


Modified 


0.54 ± 0.02 


0.46 ± 0.02 


0.51 ±0.02 



alizcd Hurst exponent of the original time series with the 
shuffled and surrogate one's, we have found that multi- 
fractality due to the broadness of the probability density 
function has more contribution than the correlation. The 
value of the Hurst exponent shows that the flow of water 
without seasonal trend is the same as random process. 
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VI. APPENDIX 



systems. We have shown that the MF-DFA1 result of the 
monthly river flows (e.g. duct water of Daugava river se- 
ries) has two crossover time scale (s x ). Our results show 
that there exists a tendency between some runoff river 
and sun activity. Indeed due to the presence of both sea- 
sonal trend and the effect of sun's period on some river 
flows, we see at least two crossover in the DFA result of 
rivers, nevertheless, it is not expected that the number 
of sunspots has superior effect on runoff water fluctua- 
tions. So this effect has different intensity on the rivers. 
These crossover time scale are due to the sinusoidal trend. 
To minimizing the effect of this trend, we have applied 
F-DFA on the river flow time series. Applying the MF- 
DFAl method on truncated data, demonstrated that the 
monthly duct water of Daugava river series is a station- 
ary time series with approximately random behavior. For 
other rivers we found long-range correlation in their sta- 
tistical behaviors. The q dependence of h(q) and r(g), 
indicated that the monthly duct water of Daugava river 
series has multifractal behavior. By comparing the gener- 



In this appendix we derive the relation between the 
exponent h(2) (DFA1 exponent) and Hurst exponent of 
a fGn signal in one dimension. We show that for such 
stationary signal the average sample variance (Eq. 4) for 
q = 2, is proportional to s h ( q \ where h(q = 2) = H. It is 
shown that the averaged sample variance F 2 {s) behaves 
as: 

= C H s 2H , (22) 
where F 2 (s\v) is defined as: 

F\s; V )= l -j2[Y v {T)-y v {i)]\ (23) 
s i=i 

and Cn is a function of Hurst exponent H . To prove 
the statement we note that the data x(k) is a fractional 



Gaussian noise (fGn), the partial sums Y(i) (Eq. 1) will 
be a ffim signal: 



= 5>(fc) - (*>• 



(24) 



fc=i 



In the DFA1. the fitting function will have the expres- 
sion (y v (i) = a v + b v i). The slope b v and intercept a v of 
a least-squares line Y(i) for every windows {y) are given 
by: 

, E'=i Y (i)i ~ i £<U Y d) E<=i * 



2 



s 3 /12 



I s b s 

a v = -V"y(i) - — Vi, 
s * — ' s * — ' 

i=l i=l 

-Eno 

c Z — ^ 



6„S 
~2~' 



(25) 



respectively. 

Using the Eqs. 4 and 25, the Eq. 23 can be written as 
follows: 

([F\s;v)]) = (±^[Y(i)-a-biA, 
^(]±Y(A + ( a i) + ^(b*) + S (ab) 



b^- 



2 ~£^) H i-X iY (i) ), (26) 



where we have discard the subscript for simplicity. The 
ffim signals is produced by using the fGn noise as follows: 



and 



Y(i) = i H x, 
Y(i)-Y(k) = \i-k\ H x, 
{{Y(z)-Y(k)] 2 ) = a 2 \i-k\ 2H , 



(27) 



(28) 



(29) 



where a 2 = (x(«) 2 ). The variance of ffim signal is: 
(Y(i) 2 ) = a 2 i 2H [30]. Expanding the left hand side of 
Eq. 29, it can be easily shown that the correlation func- 
tion of Y(i), has the following form: 



(Y{i)Y(k))=^-\i 2H + k 2H -\i-k\ 2H ]. (30) 

Finally using the Eqs. 25 and 30, it can be easily shown 
that the Eq. 26 can be written as follows: 



{[F 2 { s - v )]) v =C n ( S ) 



2H 



(31) 



To determine the C-h , we have to calculate some terms 
such as: 



i 9 ^ 

<?Y(i)Y(j)) = a - 



' j 2H + l +if H_ i \ i _ jr] 



i,3=l 



2tl 



2 s s 



i=l j=l 



tEE^ 
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a 2 / s 2//+3 g 2H+3 
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2 \2H + 2 2{2H + l) 



2 s /pi pi 

■2H+2 / / ri ^2JI^ / „./, „,\2ff. 



{l-xy H dx- / s(l-a?) a/f da; 



5+1 25+1 



(32) 



and 



E<«))4E( ! ' 2ff +i 2ff -ii-ir) 



'■7 = 1 



2W 



i=l j = l 



/ „2ff+2 f /-l \ 



a 2 s 2H+2 
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1 



^25+1 (25 + 2)(25 +i; 
Finally has the following form: 
a 2 



(33) 



(25 + 1) 25 + 2 



-3a^ 



1 



5 + 1 25 + 1 



3a 2 ^ _ 



(5 + 1) V (5+l)(25 + l) 



(34) 



Therefore the standard DFA1 exponent for a station- 
ary signal is related to its Hurst exponent as h(q = 2) = 
5. 
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